{
"cells": [
{
"cell_type": "raw",
"id": "21753ddf",
"metadata": {
"vscode": {
"languageId": "raw"
}
},
"source": [
"---\n",
"title: \"Iterative numerical methods\"\n",
"format:\n",
" revealjs:\n",
" slide-level: 3\n",
" embed-resources: true\n",
"---"
]
},
{
"cell_type": "markdown",
"id": "eeeaebbb",
"metadata": {},
"source": [
"### The nature of numerical methods"
]
},
{
"cell_type": "markdown",
"id": "e5d67be8",
"metadata": {},
"source": [
"### Example 1: Heron's method\n",
"\n",
"Objective: calculate an approximation of the square root of a number $a$\n",
"\n",
"Method:\n",
"\n",
"1. start with an initial estimate of the square root $x$\n",
"2. Correct the estimate using the formula $x_{corrected} = \\frac{1}{2} (x + \\frac{a}{x})$\n",
"3. Repeat step 2 twenty times\n",
"\n",
"The final $x$ is the square root of $a$"
]
},
{
"cell_type": "markdown",
"id": "7d2420b7",
"metadata": {},
"source": [
"### Heron's methods implementation\n",
"\n",
"Let's start with a very simple implementation..."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "8210538e",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1.414213562373095\n"
]
}
],
"source": [
"#| echo: true\n",
"\n",
"# Compute the square root of 2\n",
"a = 2.0\n",
"\n",
"# The initial estimate\n",
"x = 1.0\n",
"\n",
"# Correct 20 times, according to the formula\n",
"# x_correct = (1/2) (x + a/x)\n",
"\n",
"for i in range(20):\n",
" x = 0.5 * (x + a/x)\n",
"\n",
"print(x)"
]
},
{
"cell_type": "markdown",
"id": "479bab74",
"metadata": {},
"source": [
"### \n",
"\n",
"A few questions come to mind:\n",
"\n",
"- is the answer correct?\n",
"- why does it work? where does the \"correction\" formula come from?\n",
"- what is happening to $x$ during the corrections?\n",
"- the initial estimate $x$ matters? what if we chose a different value?\n",
"- why 20 times?\n",
"- does this \"recipe\" always work?"
]
},
{
"cell_type": "markdown",
"id": "8ee2a0c4",
"metadata": {},
"source": [
"### \n",
"\n",
"Improve the program:\n",
"\n",
"- check if the result is correct\n",
"- report what happens during the \"corrections\""
]
},
{
"cell_type": "markdown",
"id": "f16aae10",
"metadata": {},
"source": [
"### {.scrollable}"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "2048927f",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Initial estimate 1.0\n",
"Correction 1: 1.5\n",
"Correction 2: 1.4166666666666665\n",
"Correction 3: 1.4142156862745097\n",
"Correction 4: 1.4142135623746899\n",
"Correction 5: 1.414213562373095\n",
"Correction 6: 1.414213562373095\n",
"Correction 7: 1.414213562373095\n",
"Correction 8: 1.414213562373095\n",
"Correction 9: 1.414213562373095\n",
"Correction 10: 1.414213562373095\n",
"Correction 11: 1.414213562373095\n",
"Correction 12: 1.414213562373095\n",
"Correction 13: 1.414213562373095\n",
"Correction 14: 1.414213562373095\n",
"Correction 15: 1.414213562373095\n",
"Correction 16: 1.414213562373095\n",
"Correction 17: 1.414213562373095\n",
"Correction 18: 1.414213562373095\n",
"Correction 19: 1.414213562373095\n",
"Correction 20: 1.414213562373095\n",
"\n",
"Final result:\n",
"The square root of 2.0 is 1.414213562373095\n",
"The error in the result is -2.220446049250313e-16\n"
]
}
],
"source": [
"#| echo: true\n",
"#| output-location: slide\n",
"\n",
"from math import sqrt # to use with sqrt() function\n",
"\n",
"# Compute the square root of 2\n",
"a = 2.0\n",
"real_result = sqrt(a)\n",
"\n",
"# The initial estimate\n",
"x = 1.0\n",
"print(\"Initial estimate\", x)\n",
"\n",
"# Correct 20 times\n",
"for i in range(20):\n",
" x = 0.5 * (x + a/x)\n",
" print(f\"Correction {i+1}: {x}\")\n",
"\n",
"print(f\"\\nFinal result:\\nThe square root of {a} is {x}\")\n",
"print(f\"The error in the result is {x-real_result}\")"
]
},
{
"cell_type": "markdown",
"id": "29066a17",
"metadata": {},
"source": [
"### Numerical methods\n",
"\n",
"The Heron's method is a **numerical method**\n",
"\n",
"It provides a solution to a problem by *performing operations on numbers and providing a numerical solution to the problem*\n",
"\n",
"(this is increadibly vague)\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "062d5acb",
"metadata": {},
"source": [
"### Numerical methods: basic concepts\n",
"\n",
"In particular, the Heron's method works by applying a sucession of improvements called **iterations**, based on an **iteration formula**. (From now on, \"corrections\" will be called *iterations*).\n",
"\n",
"The iterations are supposed to **converge** to the solution (the method \"converges\"). After a number of sufficient iterations, the final result is an approximation to the solution of a problem.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "357402de",
"metadata": {},
"source": [
"### Different handling of convergence\n",
"\n",
"Let's modify the program to show if we can get close to the result in a number of iterations different from 20...\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "0e9abf71",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Initial estimate 1.0\n",
"Iteration 1: 1.5\n",
"Iteration 2: 1.4166666666666665\n",
"Iteration 3: 1.4142156862745097\n",
"Iteration 4: 1.4142135623746899\n",
"Iteration 5: 1.414213562373095\n",
"\n",
"Final result:\n",
"The square root of 2.0 is 1.414213562373095\n",
"The error in the result is -2.220e-16\n"
]
}
],
"source": [
"#| echo: true\n",
"#| output-location: slide\n",
"from math import sqrt # to check the answer with sqrt() function\n",
"\n",
"# Compute the square root of 2\n",
"a = 2.0\n",
"real_result = sqrt(a)\n",
"# The initial estimate\n",
"x = 1.0\n",
"print(\"Initial estimate\", x)\n",
"\n",
"tolerance = 1e-10 # max difference between x and corrected x\n",
"\n",
"# Iterations\n",
"for i in range(1000): # no more than 1000 iterations\n",
" new_x = 0.5 * (x + a/x)\n",
" print(f\"Iteration {i+1}: {new_x}\")\n",
" difference = abs(new_x - x)\n",
" x = new_x\n",
" if abs(difference) < tolerance:\n",
" break # premature exit of the for loop, before 1000 iterations\n",
"\n",
"print(f\"\\nFinal result:\\nThe square root of {a} is {x}\")\n",
"print(f\"The error in the result is {x-real_result:.3e}\")"
]
},
{
"cell_type": "markdown",
"id": "4d6971ac",
"metadata": {},
"source": [
"### The role of tolerance\n",
"\n",
"Now we had only 5 iterations instead of 20.\n",
"\n",
"The use `tolerance` is a different way to stop a numerical method that works by sucessive iterations towards a solution.\n",
"\n",
"If the \"correction amount\" is very small, the new $x$ is very close to the previous $x$ and we can consider that the approximation is sufficiently good and stop the process before reaching a large number of iterations (1000).\n",
"\n",
"`tolerance` is a small number that we set in the program and represents the maximum deviation between the new $x$ and the previous $x$ that we can consider sufficiently small to stop the iterations.\n"
]
},
{
"cell_type": "markdown",
"id": "8e9c0ddd",
"metadata": {},
"source": [
"### Keeping the \"history\" to display or plot in the end\n",
"\n",
"In a new improvement of the program, we keep the history of what is going on and display in the end only. This \"history\" will be implemented as a Python dictionary\n"
]
},
{
"cell_type": "markdown",
"id": "6ccb8e77",
"metadata": {},
"source": [
"### "
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "2019150a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Computing square root of 2.0 by Heron's method\n",
"Initial estimate: 1.0\n",
"Iteration 1: x = 1.5, dif = 0.5\n",
"Iteration 2: x = 1.4166666666666665, dif = 0.08333333333333348\n",
"Iteration 3: x = 1.4142156862745097, dif = 0.002450980392156854\n",
"Iteration 4: x = 1.4142135623746899, dif = 2.1238998197947723e-06\n",
"Iteration 5: x = 1.414213562373095, dif = 1.5949463971764999e-12\n",
"\n",
"Final result:\n",
"The square root of 2.0 is 1.414213562373095\n",
"The error in the result is -2.220e-16\n"
]
}
],
"source": [
"#| echo: true\n",
"#| output-location: slide\n",
"from math import sqrt # to check the answer with sqrt() function\n",
"\n",
"# Compute the square root of 2\n",
"a = 2.0\n",
"real_result = sqrt(a)\n",
"\n",
"# to keep the history...\n",
"history = {\"x_values\": [], \"dif_values\": [], \"initial\": 0.0, \"final\": 0.0}\n",
"\n",
"# The initial estimate\n",
"x = 1.0\n",
"history[\"x_values\"].append(x)\n",
"history[\"initial\"] = x\n",
"\n",
"tolerance = 1e-10 # max difference between x and corrected x\n",
"\n",
"# Iterations\n",
"for i in range(1000): # no more than 1000 iterations\n",
" new_x = 0.5 * (x + a/x)\n",
" difference = abs(new_x - x)\n",
" history[\"dif_values\"].append(difference)\n",
" x = new_x\n",
" history[\"x_values\"].append(x)\n",
" if abs(difference) < tolerance:\n",
" break # premature exit of the for loop, before 1000 iterations\n",
"history[\"final\"] = x\n",
"\n",
"# how many iterations\n",
"n_iter = len(history[\"dif_values\"])\n",
"\n",
"print(f\"Computing square root of {a} by Heron's method\")\n",
"print(\"Initial estimate:\", history[\"initial\"])\n",
"\n",
"\n",
"for i in range(n_iter):\n",
" iter = i+1\n",
" print(f\"Iteration {iter}: x = {history['x_values'][iter]}, dif = {history['dif_values'][i]}\")\n",
"\n",
"print(f\"\\nFinal result:\\nThe square root of {a} is {history[\"final\"]}\")\n",
"print(f\"The error in the result is {history[\"final\"] - real_result:.3e}\")"
]
},
{
"cell_type": "markdown",
"id": "6e9baaa0",
"metadata": {},
"source": [
"### Plot the history using *matplotlib*\n",
"\n",
"The history of $x$'s and $dif$'s can be plotted... \n"
]
},
{
"cell_type": "markdown",
"id": "73c395df",
"metadata": {},
"source": [
"### Preparation of a figure with two panels"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "2968cc58",
"metadata": {},
"outputs": [],
"source": [
"# don't ask...\n",
"%config InlineBackend.figure_formats = ['svg']"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "3062a9de",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#| echo: true\n",
"#| output-location: slide\n",
"from matplotlib import pyplot as plt # a necessary import to use matplotlib\n",
"\n",
"# preparation: a figure with two plots side-by-side\n",
"\n",
"f, (pleft, pright) = plt.subplots(1, 2, figsize=(8, 3.5))\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "00c94f1b",
"metadata": {},
"source": [
"### Plotting the data"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "78212659",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#| echo: true\n",
"f, (pleft, pright) = plt.subplots(1, 2, figsize=(9, 3.5))\n",
"\n",
"# plot x_values on the left, diferences on the right\n",
"\n",
"pleft.plot(range(n_iter+1), history[\"x_values\"])\n",
"pright.plot(range(1, n_iter+1), history[\"dif_values\"])\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "2b56c7dc",
"metadata": {},
"source": [
"### Improving the plots"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "0188bc7a",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#| echo: true\n",
"#| output-location: slide\n",
"# it is a bit ugly, let's improve it,\n",
"# using the power of matplotlib \"axes\" options\n",
"\n",
"# preparation\n",
"f, (pleft, pright) = plt.subplots(1, 2, figsize=(8, 3.5), tight_layout=True)\n",
"\n",
"# plot x_values on the left, differences on the right\n",
"pleft.axhline(real_result, color='darkred')\n",
"pleft.plot(range(n_iter+1), history[\"x_values\"],\n",
" marker=\"o\", markerfacecolor=\"white\")\n",
"pleft.set(xlabel=\"Iteration\", ylabel=\"$x$\",\n",
" xlim=(0, None), ylim=(1, 2), title=\"Progress\")\n",
"pleft.grid(True)\n",
"\n",
"pright.plot(range(1, n_iter+1), history[\"dif_values\"],\n",
" marker=\"o\", markerfacecolor=\"white\")\n",
"pright.set(xlabel=\"Iteration\", ylabel=\"$abs(x_{new} - x)$\",\n",
" xlim=(0, None), title=\"Corrections\")\n",
"pright.set_yscale(\"log\")\n",
"pright.grid(True)\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "e3fcc717",
"metadata": {},
"source": [
"### Exercise\n",
"\n",
"Study the effect of changing the initial estimate.\n",
"\n",
"Repeat the whole algorithm for a set of initial estimates, say seven, including negative values.\n",
"\n",
"Keep all the \"histories\", and superimpose **on the same plots** the results of the several estimations. Don't forget to add a legend for the different initial values"
]
},
{
"cell_type": "markdown",
"id": "1faec362",
"metadata": {},
"source": [
"### Solution"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "f35f6bfe",
"metadata": {},
"outputs": [],
"source": [
"#| echo: true\n",
"from math import sqrt # to check the answer\n",
"# Compute the square root of 2\n",
"a = 2.0\n",
"real_result = sqrt(a)\n",
"\n",
"# valores iniciais a testar\n",
"initial_x = [1, 3, 10, 20, -1, -3, 10, -20]\n",
"history_list = []\n",
"\n",
"for initial in initial_x:\n",
"\n",
" # to keep the history...\n",
" history = {\"x_values\": [], \"dif_values\": [], \"initial\": 0.0, \"final\": 0.0}\n",
"\n",
" # The initial estimate\n",
" x = initial\n",
" history[\"x_values\"].append(x)\n",
" history[\"initial\"] = x\n",
"\n",
" tolerance = 1e-10 # max difference between x and corrected x\n",
"\n",
" # Iterations\n",
" for i in range(1000): # no more than 1000 iterations\n",
" new_x = 0.5 * (x + a/x)\n",
" difference = abs(new_x - x)\n",
" history[\"dif_values\"].append(difference)\n",
" x = new_x\n",
" history[\"x_values\"].append(x)\n",
" if abs(difference) < tolerance:\n",
" break # premature exit of the for loop, before 1000 iterations\n",
" \n",
" history[\"final\"] = x\n",
" history_list.append(history)"
]
},
{
"cell_type": "markdown",
"id": "cf2bc97a",
"metadata": {},
"source": [
"### Solution (plots)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "10cc14fe",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#| echo: true\n",
"#| output-location: slide\n",
"# plot for all initial values\n",
"f , p = plt.subplots(1,1, figsize=(5,4))\n",
"\n",
"for history in history_list:\n",
" y = history[\"x_values\"]\n",
" x = range(len(history[\"x_values\"]))\n",
" p.plot(x, y, marker=\"o\", markersize=4,\n",
" markerfacecolor=\"white\",\n",
" label=f\"$x_0$ = {history[\"initial\"]}\")\n",
"p.set(xlabel=\"Iteration\", ylabel=\"x\", xlim=(0, None))\n",
"p.grid(True)\n",
"p.legend(ncols=2)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "04faa433",
"metadata": {},
"source": [
"### Newton's method\n",
"\n",
"- Heron's method is a particular case of the Newton's method.\n",
"\n",
"- Newton's method is a numerical method to find the zero of a function $f(x)$, from an initial estimate of the zero $x_0$ and the knowledge of the function $f(x)$ and its first derivative $f'(x)$\n",
"\n",
"- Heron's method is the Newton's method applied to function $f(x) = x^2 - a$"
]
},
{
"cell_type": "markdown",
"id": "045dd01c",
"metadata": {},
"source": [
"### Newton's method (graphical interpretation)\n",
"\n",
"As an example, how iteration $x_1$ is corrected to iteration $x_2$:\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"id": "8c46cbb5",
"metadata": {},
"source": [
"### Example 2: pH of solution of weak acid by fixed-point iteration"
]
},
{
"cell_type": "markdown",
"id": "69febd57",
"metadata": {},
"source": [
"Acetic acid: pKa = 4.76\n",
"\n",
"What is the pH of a solution of acetic acid after dissolving (and mixing) a total acid T = 0.01 M ?\n",
"\n",
"- HAc ⇄ Ac- + H+\n",
"\n",
"- $T = [Ac^-]_{eq} + [H^+]_{eq}$\n",
"\n",
"- and we will ignore the contribution of water to H+"
]
},
{
"cell_type": "markdown",
"id": "9e4a157d",
"metadata": {},
"source": [
"### \n",
"\n",
"Using fixed point iteration $x_{new} = \\sqrt{(K_a (T-x))}$"
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "2fc9d3d0",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Ka = 1.7378008287493764e-05\n",
"T = 0.01\n",
"\n",
"Initial pH estimate: 3.0\n",
"Iteration 1: pH = 3.4028787452803373\n",
"Iteration 2: pH = 3.3887621000434556\n",
"Iteration 3: pH = 3.389057710626766\n",
"\n",
"Final result:\n",
"[H+] = 4.08265e-04 M\n",
"pH = 3.39\n"
]
}
],
"source": [
"#| echo: true\n",
"#| output-location: slide\n",
"from math import log10\n",
"\n",
"pKa = 4.76\n",
"T = 0.01 # M (dissolved)\n",
"\n",
"# Compute Ka\n",
"Ka = 10**(-pKa)\n",
"print(\"Ka = \", Ka)\n",
"print(\"T = \", T)\n",
"\n",
"# x, the [H+] = [Ac-], verifies Ka = x**2 / (T - x)\n",
"# or\n",
"# x = sqrt(Ka * (T - x))\n",
"\n",
"# to keep the history...\n",
"history = {\"x_values\": [], \"pH_values\": [], \"deltapH\": []}\n",
"\n",
"\n",
"# The initial estimate\n",
"x = 0.001\n",
"history[\"x_values\"].append(x)\n",
"history[\"pH_values\"].append(-log10(x))\n",
"\n",
"tolerance = 0.01 # max difference between pH and corrected pH\n",
"\n",
"# Iterations\n",
"for i in range(100): # no more than 100 iterations\n",
" old_pH = -log10(x)\n",
"\n",
" new_x = (Ka * (T-x))**0.5 # the correction\n",
"\n",
" pH = -log10(new_x)\n",
" history[\"pH_values\"].append(pH)\n",
"\n",
" dif_pH = abs(pH - old_pH)\n",
" x = new_x\n",
" history[\"x_values\"].append(x)\n",
" if abs(dif_pH) < tolerance:\n",
" break # premature exit of the for loop\n",
"\n",
"# how many iterations\n",
"n_iter = len(history[\"x_values\"]) -1\n",
"\n",
"print(\"\\nInitial pH estimate:\", history[\"pH_values\"][0])\n",
"\n",
"for i in range(n_iter):\n",
" iter = i+1\n",
" print(f\"Iteration {iter}: pH = {history['pH_values'][iter]}\")\n",
"\n",
"print(f\"\\nFinal result:\\n[H+] = {history[\"x_values\"][-1]:.5e} M\")\n",
"print(f\"pH = {history[\"pH_values\"][-1]:.2f}\")\n"
]
},
{
"cell_type": "markdown",
"id": "2494da08",
"metadata": {},
"source": [
"### \n",
"\n",
"Notice how tolerance now applies to pH values and not to $x$, the concentration of H+ = concentration of Ac-"
]
},
{
"cell_type": "markdown",
"id": "55ffc515",
"metadata": {},
"source": [
"### Exercise\n",
"\n",
"Plot the pH as a function of dissolved HAc,\n",
"\n",
"in the range 0.0001 M to 0.1 M"
]
},
{
"cell_type": "markdown",
"id": "92cedadd",
"metadata": {},
"source": [
"### Solution\n",
"\n",
"(of exercise, not of acetic acid)"
]
},
{
"cell_type": "code",
"execution_count": 55,
"id": "97c0a1e2",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" T pH\n",
"0.00010000 4.47\n",
"0.00014384 4.38\n",
"0.00020691 4.28\n",
"0.00029764 4.20\n",
"0.00042813 4.11\n",
"0.00061585 4.02\n",
"0.00088587 3.94\n",
"0.00127427 3.85\n",
"0.00183298 3.77\n",
"0.00263665 3.69\n",
"0.00379269 3.61\n",
"0.00545559 3.52\n",
"0.00784760 3.44\n",
"0.01128838 3.36\n",
"0.01623777 3.28\n",
"0.02335721 3.20\n",
"0.03359818 3.12\n",
"0.04832930 3.04\n",
"0.06951928 2.96\n",
"0.10000000 2.88\n",
"\n",
"Is equilibrium equation verified?\n",
"(x*x / (T-x)) / Ka = 0.999\n",
"(x*x / (T-x)) / Ka = 1.001\n",
"(x*x / (T-x)) / Ka = 1.000\n",
"(x*x / (T-x)) / Ka = 1.000\n",
"(x*x / (T-x)) / Ka = 1.000\n",
"(x*x / (T-x)) / Ka = 1.000\n",
"(x*x / (T-x)) / Ka = 1.000\n",
"(x*x / (T-x)) / Ka = 1.000\n",
"(x*x / (T-x)) / Ka = 1.000\n",
"(x*x / (T-x)) / Ka = 1.000\n",
"(x*x / (T-x)) / Ka = 1.000\n",
"(x*x / (T-x)) / Ka = 1.000\n",
"(x*x / (T-x)) / Ka = 1.000\n",
"(x*x / (T-x)) / Ka = 1.000\n",
"(x*x / (T-x)) / Ka = 1.000\n",
"(x*x / (T-x)) / Ka = 1.000\n",
"(x*x / (T-x)) / Ka = 1.000\n",
"(x*x / (T-x)) / Ka = 1.000\n",
"(x*x / (T-x)) / Ka = 1.000\n",
"(x*x / (T-x)) / Ka = 1.000\n"
]
}
],
"source": [
"#| echo: true\n",
"#| output-location: slide\n",
"from math import log10\n",
"import numpy as np\n",
"\n",
"pKa = 4.76\n",
"Ka = 10**(-pKa)\n",
"tolerance = 0.001 # max difference between pH and corrected pH\n",
"\n",
"# now T varies\n",
"T_values = np.logspace(-4, -1, 20) # 20 points, from 1e-4 to 0.1\n",
"pH_values = [] # to hold the computed values\n",
"\n",
"for T in T_values:\n",
" # no need to keep the history, only the pH\n",
" # The initial estimate\n",
" x = T/100\n",
"\n",
" # Iterations\n",
" for i in range(100): # no more than 100 iterations\n",
" old_pH = -log10(x)\n",
" new_x = (Ka * (T-x))**0.5 # the correction\n",
" pH = -log10(new_x)\n",
"\n",
" dif_pH = abs(pH - old_pH)\n",
" x = new_x\n",
" if abs(dif_pH) < tolerance:\n",
" break # premature exit of the for i loop\n",
" pH_values.append(pH)\n",
"\n",
"print(\" T pH\")\n",
"for T, pH in zip(T_values, pH_values):\n",
" print(f\"{T:10.8f} {pH:.2f}\")\n",
"\n",
"print(\"\\nIs equilibrium equation verified?\")\n",
"for T, pH in zip(T_values, pH_values):\n",
" x = 10**(-pH)\n",
" eq_ratio = (x*x / (T-x)) / Ka\n",
" print(f\"(x*x / (T-x)) / Ka = {eq_ratio:.3f}\")"
]
},
{
"cell_type": "markdown",
"id": "0da0aeff",
"metadata": {},
"source": [
"### Plots"
]
},
{
"cell_type": "code",
"execution_count": 54,
"id": "8239c3bf",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"f, ax = plt.subplots(1,1, figsize=(6,4))\n",
"ax.plot(T_values, pH_values, color=\"darkred\")\n",
"ax.set(xlabel=\"Molarity\", ylabel=\"pH\")\n",
"ax.set_xscale(\"log\")\n",
"ax.grid()\n",
"plt.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "PMN",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.12"
}
},
"nbformat": 4,
"nbformat_minor": 5
}